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1  Introduction 


This  third  semi-annual  progress  report  contains  a  summary  of  work  accomplished  on  O. 
N.  R.  contract  number  N00014-86-K.-0370,  Delay-Doppler  Radar-Imaging ,  during  the  period 
from  1  June  1987  to  30  November  1987.  This  six-month  period  is  the  last  of  three  making  up 
the  18  month  period  of  this  contract.  A  one-year  continuation  contract  entitled  High  Resolution 

Radar-Imaging  has  been  awarded  for  the  period  1  December  1987  to  30  November  1988. 

/ 

The  goal  of  this  project  is  to  formulate  and  investigate  new  approaches  for  forming 
images  of  radar  targets  from  spotlight-mode,  delay-doppler  measurements.  These  measurements 
could  be  acquired  with  a  high-resolution  radar-imaging  system  operating  with  an  optical-  or 
radio-frequency  carrier.  Two  approaches  are  under  study.  The  first  is  motivated  by  an 
image-reconstruction  algorithm  used  in  radionuclide  imaging  called  the  confidence-weighted 
algorithm ;  here,  we  will  refer  to  this  approach  as  the  chirp-rate  modulation  approach.  The 
second  approach  is  based  on  more  fundamental  principles  which  starts  with  a  mathematical 
model  that  accurately  describes  the  physics  of  an  imaging  radar-system  and  then  uses 
statistical-estimation  theory  with  this  model  to  derive  processing  algorithms;  we  will  refer  to 
this  as  the  estimation-theory  approach. 

Spotlight-mode  high-resolution  radar-imaging  relies  upon  the  relative  motion  between  the 
transmitter,  target,  and  receiver.  The  target  is  illuminated  by  a  series  of  transmitted  pulses. 
The  return  for  each  pulse  is  a  superposition  of  reflections  from  various  locations  on  the  target, 
with  each  location  affecting  the  pulse  by  introducing  a  time  delay,  doppler  shift,  and 
reflectance  gain.  The  returns  are  processed  to  produce  an  image  of  the  target.  Two  types  of 
images  are  possible.  One  is  a  map  in  delay/doppler  or  range/cross-range  coordinates  of  the 
target’s  complex-valued  reflectivity,  which  indicates  the  amplitude-gain  applied  to  the  incident 
radar-pulse  by  each  location  on  the  target.  The  other  is  a  map  in  the  same  coordinates  of  the 
target’s  scattering  function,  which  indicates  the  power-gain  applied  at  each  location  on  the 

target. ,  We  are  studying  approaches  for  producing  both  types  of  images. 

/ 

A  common  approach  discussed  in  the  literature  and  implemented  in  practice  is  to  use  the 
same  transmitted  pulse  for  each  illumination  of  the  target.  The  returns  are  processed  using  a 
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two-dimensional  Fourier  transform  to  produce  the  target’s  image  in  delay/doppler  or 
range/cross-range  coordinates  [1,2].  One  of  our  goals  is  to  compare  images  produced  in  this 
standard  way  with  those  produced  using  the  alternative  approaches  we  are  developing. 

Bernfeld  [3]  appears  to  be  the  first  to  introduce  the  idea  for  radar  imaging  of  modifying 
the  pulse  shape  on  successive  illuminations  of  the  target.  We  are  using  Bernfeld’s  idea  in  the 
chirp-rate  modulation  approach.  With  this,  the  FM  chirp-rate  of  each  pulse  is  varied  so  that 
the  angles  made  by  the  ambiguity  functions  in  the  delay/doppler  plane  are  caused  to  vary  over 
the  full  range  of  angles  between  0°  and  180°.  Use  is  then  made  of  the  fact  that,  on  the  average, 
the  square-law  envelope-detected  output  of  a  receiver  consisting  of  a  bank  of  filters  matched  to 
doppler  shifted  versions  of  a  transmitted  pulse  is  the  two-dimensional  convolution  of  the 
ambiguity  function  of  the  pulse  with  the  scattering  function  of  the  target  [4],  an  output  we  call 
the  delay/doppler  power  function.  Given  the  delay/doppler  power  functions  for  each 
illumination,  the  target’s  scattering  function  can  be  determined  using  the  confidence- weighted 
algorithm  [4]. 

Statistical  models  for  radar  echoes  from  spatially  extended  targets  that  are  rough  compared 
to  a  wavelength  of  the  radar  carrier  are  given  by  Van  Trees  [5]  for  microwave  frequencies  and 
Shapiro,  Capron,  and  Harney  [6]  for  optical  frequencies.  In  the  estimation-theory  approach,  we 
are  using  these  models  with  statistical  estimation-theory  to  develop  new  algorithms  for 
producing  high  resolution  images  of  radar  targets  from  delay/doppler  data  [7], 

Work  accomplished  during  the  reporting  period  is  summarized  in  the  following  section. 


2  Summary  of  Work  Accomplished 

2.1  Estimation-Theory  Approach  to  Imaging 

During  this  reporting  period,  reference  [7,  see  Appendix  1  for  a  preprint]  was  revised  and 
submitted  to  the  IEEE  Transactions  on  Information  Theory  for  review  and  publication.  This 
manuscript  contains  a  discussion  of  the  mathematical  model  we  are  using  to  describe  radar 
returns  from  diffuse  radar  targets.  Equations  are  derived  for  the  maximum-likelihood  estimate 
of  the  target’s  reflectivity  process  and  its  scattering  function. 


In  this  project,  the  complex  envelope  of  the  radar-echo  data  is  described  by: 

r  ( t )  -  \j2Fj  j '  s(f-r)b^t-^r.rjdr  +  iti(t). 

where  £t  is  the  transmitted  energy,  s(t)  is  the  complex  envelope  of  the  transmitted 
pulse-sequence,  w{t)  is  broad  band  Gaussian  noise,  and  b(tj)  is  the  target’s  reflectivity  at  time  t 
for  all  the  reflecting  patches  on  the  target  at  two-way  delay  r.  This  model  describes  the  radar 
data  as  a  superposition  of  echoes  from  all  reflecting  patches  on  the  target  plus  an  additive  noise. 
The  reflectivity  b(tj)  is  assumed  to  be  a  zero-mean,  wide-sense  stationary,  uncorrelated-scat- 
terers  (WSSUS)  Gaussian  process  with  a  scattering  function  S(/,t)  in  doppler  /  and  delay  r 
coordinates;  this  model  for  the  reflectivity  is  motivated  in  references  [5]  and  [6].  The  scattering 
function  is  the  power-density  spectrum  of  the  reflectivity  process  b(tj)  at  two-way  delay  r. 
The  problem  we  address  in  [7]  is  that  of  forming  the  maximum-likelihood  estimate  of  the 
scattering  function  based  on  data  r(t).  This  problem  is  addressed  in  a  general  framework,  for 
example  with  the  form  of  the  transmitted  pulse  left  arbitrary,  so  that  the  fundamental  equations 
describing  the  estimate  can  be  obtained  for  a  wide  variety  of  potentially  interesting  special 
cases.  The  estimate  satisfies  a  nonlinear  integral-equation;  we  describe  an  iterative  approach  for 
solving  the  equation  numerically.  The  unique  feature  of  our  approach  is  that  we  have  addressed 
the  problem  of  radar  imaging  by  starting  with  a  basic  statistical  model  for  the  return  data  and 
then  deriving  an  image-formation  algorithm  using  estimation  theory.  We  expect  that  improved 
target  images  will  be  obtained  in  those  situations  where  the  reflectivity  model  is  an  accurate 
description  of  reality. 

At  the  present  time,  we  are  beginning  to  investigate  scattering-function  estimates  in 
computer  simulations,  and  we  are  attempting  to  extend  the  model  of  the  target's  reflectivity  to 
include  specular  or  glint  components  in  the  echo  signal  in  addition  to  the  diffuse  component 
contained  in  the  model  of  [7],  We  also  wish  to  incorporate  constraints,  implied  by  targets  of 
interest,  into  the  image-formation  process. 

2.2  Chirp-Rate  Modulation  Approach  to  Imaging 

The  focus  of  our  research  on  the  chirp-rate  modulation  approach  during  the  reporting 
period  has  been  the  analytical  investigation  of  a  linear  architecture  based  on  a  bank  of  bandpass 


matched-filters.  The  motivation  for  this  resulted  from  a  concern  that  the  nonlinear  architecture 
being  used  in  our  investigations,  a  bank  of  bandpass  matched-filters  followed  by  a  square-law 
envelope  detector,  results  in  a  loss  of  coherence  in  the  processing  of  data  from  successive  target 
illuminations  and  that  this  loss  could  result  in  degraded  target  images. 

We  consider  the  above  model  for  return  data  r{t),  but  temporarily  without  the  additive 
noise  w(t).  By  writing  the  reflectivity  process  b(tj)  in  terms  of  its  Fourier  transform  c(fdS), 

b(t.T)-  f  c(/d , r )e‘2nl df  d  , 


we  can  express  the  return  data  according  to 


r(0=  c{f  d,T  )s(t  -  r  )e  df  ddx 


The  bank  of  bandpass  matched  filters  produces  the  complex-valued  delay/doppler  image  q(t,f) 
defined  by 


<7(1 .  /)  *  r(u)s*(u-  t)e 


■  l2n  /  (u- I) 


c(/d.r)<^T-f,/d-/je 


df  dd  re' 


where 


t  .  \  I  1  \  [  \  \  ,2  n,dz  , 

♦  ,--r  2*-T  e  * 


is  the  complex- valued  time-frequency  autocorrelation  function  of  the  transmitted  pulse  (also 
sometimes  called  the  complex-valued  ambiguity  function).  Thus,  the  data  at  the  output  of  the 
bank  of  matched  filters  are  a  linear  functional  of  c(/d,r)  and  the  time-frequency  autocorrelation 
function.  The  goal  is  to  solve  the  inverse  problem  of  determining  c(f&,T)  from  a  collection  of 
filter  outputs  qfrj),  where  the  /th  output  results  from  the  <th  transmitted  pulse  in  the 
chirp-rate  modulated  sequence.  One  potentially  interesting  observation  that  can  be  made  is  that 
if  the  time-frequency  autocorrelation  function  is  slowly  varving  (ideally,  a  constant)  over  the 
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extent  of  the  target  image  c(/d,r)  in  the  delay/doppler  plane,  then  <j(/,r)  is  related  to  the 
two-dimensional  Fourier  transform  of  c(/d,r),  and  the  reflectivity  image  may  then  be  formed  by 
inverse  Fourier-transformation  in  the  absence  of  noise. 

In  the  progress  report  submitted  for  the  first  six  month  period  of  this  contract,  there  was 
a  table  listing  various  parameter  values  being  used  in  computer  simulation  experiments  with  the 
confidence-weighted  algorithm.  An  update  of  this  table  is  contained  in  Appendix  2. 
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ABSTRACT 

* 

In  this  paper,  we  present  a  new  approach  to  high-resolution  radar  imaging.  The  starting 
point  is  a  model  of  the  radar  echo-signal  based  on  the  physics  governing  radar  reflections.  This 
model  has  been  used  several  times  in  the  past  for  describing  radar  targets  that  are  rough 
compared  to  the  wavelength  of  the  transmitted  radiation.  Without  specifying  precisely  what  the 
transmitted  signal  is,  we  derive  a  general  estimation-based  procedure  for  obtaining  images. 
After  discretizing  the  model,  the  radar  imaging  problem  reduces  to  the  task  of  estimating 
discretized  second-order  statistics  of  the  reflectance  process  of  the  target.  A  maximum-likeli¬ 
hood  estimate  of  these  statistics  is  obtained  as  the  limit  point  of  an  expectation-maximization 
algorithm. 


t  This  work  was  supported  by  the  Office  of  Naval  Research  under  contract  N00014-86-K0370. 


1.  Introduction 


Radar  systems  can  be  used  to  produce  high-resolution  images  of  a  reflecting  target.  This  is 
accomplished  by  illuminating  the  target  with  a  series  of  pulses  and  observing  the  return  echos. 
Each  patch  on  the  target  introduces  a  certain  amount  of  propagation  delay  and  doppler  shift  to 
a  pulse  it  reflects,  me  amount  depending  on  the  range  and  velocity  of  the  patch  relative  to  the 
radar  transmitter  and  receiver.  The  beam-width  of  the  radar  antenna  relative  to  the  si2e  of  the 
target  is  an  important  consideration.  Images  can  be  produced  by  scanning  a  narrowly  focused 
beam  over  the  target  in  some  type  of  raster  pattern  and  then  displaying  the  received  power  in 
delay  and  doppler  or,  equivalently,  range  and  cross-range  coordinates.  Images  can  also  be 
formed  by  illuminating  the  entire  target  in  spotlight  mode  with  a  wide,  relatively  unfocused 
beam.  The  received  signal  for  each  illumination  is  then  a  complicated  superposition  of  the 
echos  received  from  all  the  patches  that  make  up  the  extended  surface  of  the  radar  target.  Our 
concern  will  be  with  forming  images  of  rotating,  rough  targets  using  a  spotlight-mode  radar. 

We  will  denote  the  complex  envelope  of  the  signal  transmitted  by  the  radar  as 
(2£t)1/2St(/),  where  E t  is  the  transmitted  energy,  and  st(i)  is  normalized  to  unit  energy.  The 
particular  form  of  this  signal  will  not  need  to  be  specified.  The  expressions  we  obtain  for 
producing  an  image  can  then  be  specialized  for  any  signal  of  choice,  including  the 
stepped-frequency  and  wideband  chirp  waveforms  used  in  practice,  as  discussed  by  Wehner  [1] 
and  Mensa  (2]  and  the  chirp-rate  modulated  waveforms  discussed  by  Bernfeld  [3]  and  Snyder, 
Whitehouse,  Wohlschlaeger  and  Lewis  [4],  When  specializing  .st(').  it  should  be  kept  in  mind 
that  this  represents  the  entire  sequence  of  transmitted  pulses  that  illuminate  the  target. 

Following  Walker  [5],  consider  a 

small,  nonfluctuating  reflector  that  is  * 

rotating  about  an  axis  at  the  rate  of  fr  -  • 

revolutions  each  second,  as  shown  in  Fig.  R  ,  ir  ;f 

1.  The  distance  from  the  radar  tnasmit-  radar 

T/R 

ter/receiver  to  the  axis  of  rotation  is  R o, 
and  the  distance  to  the  reflector  at  time  t 

is  given  approximately  by  Figure  1.  Radar  geometry  for  a  small  rotating 

reflector 


/?(f)  =  R  o  +  x0cos(27r  /,f)  +  y0sin(2/r/r/!- 

provided  flo  »  (.vo2  +  yo2)1/2.  Then,  the  radar  echo-signal  sr(0  received  following  an 
illumination  by  $x(0  WU1  be  of  the  form 

s*(0  =  /2FrsT(!-r)6 

where  r  =  2R(t)/c  is  the  two-way  propagation  delay  to  the  reflector,  with  c  being  the 
propagation  velocity.  The  quantity  b  is  a  complex-valued  scale  factor  which  models  the 
strength  of  the  received  echo.  This  scale  factor  will  be  called  the  reflectivity.  It  includes  the 
effects  of  inverse  square-law  attenuation  experienced  by  the  propagating  radiation  and, 
importantly,  the  properties  of  the  reflector  that  are  significant  in  the  electromagnetic  scattering 
interaction,  including  its  shape,  size  and  surface  properties.  More  generally,  the  reflectivity  can 
vary  with  time  because  the  aspect  and,  therefore,  the  scattering  interaction  with  the  reflector 
will  vary  as  it  rotates.  If,  as  discussed  by  Walker  [5],  the  radar  data  5r(/)  are  examined  over  a 
small  interval  of  time,  then  the  delay  r  and  doppler  shift  /d  can  be  approximated  by 

t  =  2c~'(/?0+.y0;, 

and 


r  2 clR  .2  0  , 


where  X  is  the  wavelength  at  the  carrier  frequency  of  the  radar.  Thus,  the  reflectivity  b ,  range 
x0,  cross-range  .y0,  relative  to  the  coordinate  axis,  can  be  determined  from  the  amplitude,  delay, 
and  doppler  information  contained  in  the  radar  data.  Extracting  this  information  permits  an 
image  of  the  reflector  to  be  formed  by  displaying  |  b  |  or  |  b  |2  at  the  appropriate  location  in 
range  and  cross-range  coordinates.  The  maximum  delay  and  doppler  shift  are  determined  by 


the  distance  (.v02  +  yo2)1!2  of  the  reflector  from  the  coordinate  center  about  which  it  rotates  and 
the  rotation  rate  /r  ;  more  generally,  the  extent  of  a  reflector  in  delay  and  doppler  is  determined 
by  the  physical  extent  of  the  reflector  and  the  rotation  rate. 


Now  consider  a  spatially  extended  target  that  is  rotating.  A  patch  on  the  surface  with  a 
two-way  delay  in  the  interval  (r,r+Ar)  reflects  a  signal  that  is  incident  on  the  patch  at  time  t 
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with  a  reflectance  strength  biu)Ar.  Consequently,  the  complex  envelope  of  the  received 
echo-signal  5r(/)  following  the  illumination  of  the  target  by  s^it)  is  given  by  the  following 
superposition  of  returns  from  reflecting  patches  at  all  the  two-way  delays  r  : 


s#(0  =  \2Ft  j"  sT(t-r)b^t-^r,rjdr.  (1) 

-  ao 

The  total  received  signal,  r(t),  is  also  assumed  to  be  corrupted  by  an  independent,  additive 
noise, 

r(0“  sR(t)*w(kt).  (2) 

where  w(/)  is  a  complex-valued,  white  Gaussian-process  with  a  mean  of  zero  and  a  covariance 
function  defined  by 

E [  tt  ( t ) u,  * (  f  )  ]  =  X  06(t  -  t '  )  ,  (3) 

where  "  *  "  denotes  complex  conjugation.  We  refer  to  b(tj)  as  the  reflectance  process.  This  is  a 
complex-valued  random  process. 

There  are  two  images  which  may  be  displayed  as  the  result  of  processing  /•(/)  with  an 
estimation  procedure.  One  is  an  estimate  of  the  reflectance  process  b(tj)  itself  and  the  other  is 
an  estimate  of  the  covariance  or,  equivalently,  the  spectral  density  of  this  process.  These  may 
be  regarded  as  conditional  first-order  and  a  second-order  statistics  of  the  reflectivity  process, 
respectively,  in  terms  of  the  radar  data  (2). 

If  b(fj)  is  deterministic,  define  c(r,r)  to  be  its  Fourier  transform  in  the  t  variable, 

c(  /  .  t  )  -  f  b(t  .r)<?" t2n"  ell .  (4) 

-  w 

The  function  c(/,r)  then  contains  information  about  the  target  in  delay,  r,  and  doppler,  f, 
coordinates.  An  image  of  the  target  is  obtained  by  placing  the  magnitude  or  squared  magnitude 
of  this  function  into  delay  and  doppler  bins,  We  refer  to  this  as  the  reflectance  image.  This 
transform  can  be  obtained  in  a  variety  of  ways,  depending  on  the  signal  sj(t)  selected  to 
illuminate  the  target.  For  the  stepped-frequencv  signals  used  in  practice,  the  usual  approach 
consists  of  two  operations  described  by  Wehner  [1],  The  first  is  to  place  the  data  into  delay  (or, 
range)  bins  bv  separately  Fourier  transforming  V  sample  values  of  the  received  signal  acquired 


for  each  transmitted  group  of  /V  stepped-frequency  pulses.  The  resulting  delay-binned  data  are 
placed  in  the  rows  of  an  NxN  matrix,  where  each  row  contains  the  transformed  data  from  one 
pulse-group.  In  the  second  operation,  the  columns  of  this  matrix  are  Fourier  transformed  to 
obtain  a  doppler  (or,  cross-range)  profile  at  each  delay.  The  resulting  two-dimensional  array  is 
intended  to  be  a  discrete  version  of  c(/,r)  in  delay  (range)  and  doppler  (cross-range)  coordinates. 
This  processing  based  on  two-dimensional  Fourier  transforms  is  derived  using  a  strictly 
deterministic  analysis  and  so  does  not  account  for  statistical  properties  of  the  reflectance  process 
or  for  noise  that  may  be  present.  A  similar  processing  is  employed  for  the  linear  FM-chirp 
signals  also  used  in  practice  for  radar  imaging  [1,  2]. 

For  situations  in  which  the  target’s  surface  is  rough  compared  to  the  wavelength  at  the 
carrier  frequency,  b(tj)  may  be  taken  to  be  a  complex- valued  Gaussian  random  process,  as 
discussed  by  Shapiro,  Capron,  and  Harney  [6]  for  radar  systems  operating  at  laser  frequencies 
and  Van  Trees  [7]  at  microwave  frequencies.  If  there  are  no  glint  or  specular  components  in 
the  echo,  then  this  is  a  zero  mean  process  with  covariance 


E[6(t,T)b*(I',r/)]  =  A'(I-I',r)6(T-r'). 


The  delta  function  in  this  expression  results  from  postulating  that  each  reflecting  patch 
introduces  an  uncorrelated  contribution  to  the  echo  signal.  That  the  function  K(t-t  *,r)  depends 
only  on  the  difference  of  i  and  and  not  on  t  and  l'  separately,  results  from  postulating  that 
the  reflectance  process  is  wide-sense  stationary  for  each  delay.  A  reflectance  process  with  these 
properties  is  said  by  Van  Trees  [7]  to  possess  wide-sense  stationary,  uncorrelated  scatterers 
(WSSUS).  Assuming  that  the  reflectance  process  has  these  properties,  the  delay-doppler  data 
associated  with  the  radar  target  may  be  obtained  from  the  Fourier  transform  of  K(tj)  in  the  t 
variable. 


S(/-r)“J  A'(/ ,r)exp(- j2n  ft)dl .  (6) 

-  oo 

The  function  Sifj)  is  called  the  scattering  function  of  the  target  and,  as  a  function  of  /,  is  the 
power-density  spectrum  of  the  reflectance  process  at  each  delay  r.  S(/.r)A/Ar  is  the 
mean-square  strength  (or  power)  of  the  reflectance  of  all  patches  on  the  target  having  a  doppler 
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shift  in  the  interval  [/,/VA/)  and  a  delay  in  the  interval  [r,r+Ar).  The  scattering  function  may  be 
viewed  in  delay  and  doppler  coordinates  as  an  image  of  the  target.  We  call  this  the  scattering 
function  image. 

Our  approach  to  forming  radar  images  will  be  to  use  maximum-likelihood  methods  with 
the  data  model  in  (2)  to  estimate  the  scattering  function.  We  will  also  obtain  an  estimate  of  the 
reflectance  process.  Model  based  approaches  that  use  statistical  estimation-theory  techniques  to 
derive  image-formation  algorithms  appear  less  frequently  in  the  large  literature  about  radar 
imaging  than  do  the  deterministic  approaches  outlined  above.  One  example  is  that  of  Frost, 
Stiles,  Shanmugan,  and  Holtzman  [8],  who  use  a  multiplicative  model  and  Wiener-filtering 
techniques.  The  approach  we  will  describe  below  differs  in  that  the  model  (2)  we  adopt  of  the 
echo  signal  is  more  complicated  than  a  simple  multiplicative  one  and  depends  explicitly  on  the 
transmitted  waveform  through  a  spatial  integration  over  the  reflecting  target.  We  also  do  not 
restrict  the  processing  to  be  linear;  in  particular,  we  show  that  the  processing  of  the  received 
data  for  producing  the  maximum-likelihood  estimate  of  the  scattering  function  and  a 
corresponding  estimate  of  the  reflectance  process  is  not  linear.  An  approach  for  estimating 
scattering  functions  of  spread  channels  is  given  by  Gaarder  [9],  who  cites  earlier  work  on  the 
subject  by  Green  [10],  Kailath  [11],  Gallager  [12],  Hagfors  [13,14],  Price  [15],  Levin  [16], 
Abraham  [17],  Sifford  [18],  and  Reiffen  [19].  Gaarder  assumes  a  specific  processing 
architecture  in  the  form  of  a  cascade  of  a  linear  filter,  square-law  envelope  detector,  and 
another  linear  filter,  and  he  claims  that  this  processing  is  either  more  general  than  or  equivalent 
to  those  of  most  previous  authors.  Our  approach  differs  in  that  no  particular  processing  is 
assumed  in  advance;  rather,  we  derive  the  processing  to  produce  the  estimates,  starting  from  a 
model  for  the  received  data  and  applying  recent  results  in  maximum-likelihood  estimation.  The 
processing  which  results  is  quite  distinct  from  that  discussed  by  Gaardner. 

For  our  new  approach  to  radar  imaging,  we  adopt  the  WSSUS  model  of  a  diffuse  radar 
target  described  by  Shapiro,  Capron,  and  Harney  [6]  and  Van  Trees  [7],  We  treat  both  the 
reflectance  process  and  its  second-order  statistic,  the  scattering  function,  as  unknown  quantities. 
The  iterative  approach  we  develop  for  forming  images  yields  the  maximum-likelihood  estimate 
of  the  scattering  function  and,  simultaneously,  the  conditional-mean  estimate  of  the  reflectance 
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process  based  on  statistics  which  are  consistent  with  the  estimated  scattering  function.  Thus, 
both  of  the  quantities  treated  separately  in  other  imaging  schemes  are  produced  simultaneously 
with  our  new  approach,  which  is  a  unique  and  important  aspect  of  our  approach. 


We  will  develop  a  necessary  condition,  called  the  trace  condition ,  which  the  maximum-like¬ 
lihood  estimate  of  the  target’s  scattering  function  must  satisfy.  This  equation  appears  to  be  very 
hard  to  solve  analytically.  As  a  consequence,  we  reformulate  the  imaging  problem  using  the 
concept  of  incomplete-complete  data  spaces  and  then  use  the  expectation-maximization 
algorithm  of  Dempster,  Laird,  and  Rubin  [20]  to  derive  an  iterative  algorithm  for  producing  the 
maximum-likelihood  estimate  of  the  scattering  function.  The  technique  we  use  to  accomplish 
this  parallels  that  described  by  Miller  and  Snyder  [21]  for  power-spectrum  estimation  and 
extends  their  work  to  include  indirect  measurements  of  the  process  whose  spectrum  is  sought; 
the  process  is  now  measured  following  the  linear  transformation  and  additive  noise  seen  in  (2). 
As  shown  by  Turmon  and  Miller  [22],  this  is  a  high-resolution  approach  to  spectrum  estimation 
which  results  in  estimated  spectra  with  less  bias  and  mean-square  error  than  other  recently 
developed  approaches  discussed  in  the  literature.  We  expect  that  similar  improvements  will  be 
seen  in  radar  images  of  scintillating,  diffuse  targets  when  this  new  technique  is  used. 


2.  Maximum- Likelihood  Imaging  for  the  Incomplete-Data  Model 


For  reasons  that  will  become  evident  in  the  next  section,  we  term  the  data  r(t)  in  (2)  the 
incomplete  data  for  the  radar-imaging  problem.  The  model  given  in  the  Introduction  for  these 
data  is  as  the  sum  of  the  radar  echo-signal  sr(/)  of  (1)  and  an  independent,  white  noise-process 
w(t).  We  may,  therefore,  state  the  problem  of  imaging  a  diffuse  radar-target  as  that  of 
estimating  the  scattering  function  5(/.r)  or,  equivalently,  the  covariance  function  K(t.r)  given 
radar-return  data  {r(/),  Tx  <  t  <  Tf)  on  an  observation  interval  {Tx,Tf).  In  this  section,  we  first 
discretize  the  model  for  the  incomplete  data,  and  then  we  derive  and  discuss  a  necessary 
condition,  called  the  trace  condition ,  which  the  maximum-likelihood  estimate  of  the  discretized 
version  of  K(t.r)  must  satisfy. 


discrete  model 

In  anticipation  of  using  discrete-time  processing  of  radar  data  to  produce  images,  we  now 
state  the  discrete  version  of  our  model  as  follows.  We  are  given  N  samples  of  the 
complex-valued  radar-data  corresponding  to  (2), 

r(n)  =  sR(n)  +  w(n),  n  =  0, 1 . V  -  1  ,  (7) 

where  w(n)  is  a  white  Gaussian-sequence  with  zero  mean  and  covariance 

E[iu(n)w*  (n' )]  =  A'  0  6  „  „  ,  (8) 

and  where  the  signal  samples  corresponding  to  (1)  are  given  by 

*  30 

sR(n)  =  J2ET  Y,  sT(n.t)b(n,i).  n- 0.1 . V  -  1  .  ((?) 

In  this  expression,  we  define  sj(n.i)  and  b{n.t)  in  terms  of  the  transmitted  signal  and  the 
reflectance  process,  respectively,  according  to: 

sT(n.i)  =  sT(n.M-i  Jr). 
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where  A t  and  A r  are  the  sampling  intervals  adopted  in  the  discretization  in  time  and  delay, 
respectively.  We  assume  that  the  target  has  a  finite  extent  in  range ,  so  bOu)  and,  therefore 
also,  terms  forming  the  sum  in  (9)  are  zero  for  i  outside  the  /r  (here,  the  subscript  R  denotes 

range )  values  m.  m+1 .  m+I^-  1  starting  from  the  minimum  two-way  delay  corresponding  to  m. 

This  discrete  reflectance  is  a  Gaussian  sequence  with  zero  mean  and  covariance  given  by 

E[b(n,i)b*(n',t')]  =  K(n-n',i)6tl.,  <■-) 

where  Sy  is  the  Kronecker  delta-function.  The  discrete  scattering  function  S(f.i)  is  the  Fourier 
transform  of  K(n.i)  in  the  n  variable. 


S (  v  ,  ( )  =  V  A' (n  ,  t ) e x p ( -  j2nv n) . 

n-  -  ® 


(13) 


The  imaging  problem  for  the  discrete  model  is  to  estimate  S( ir  .;),  or  equivalently  the  covariance 
function  K(n,i ),  for  all  frequencies  u  spanning  the  target  in  doppler,  and  for  all  delays  i 
spanning  the  target  in  delay,  given  the  radar  data  { r(n ),  n=0.  1 . V-/). 


matrix  model 

These  discrete  equations  may  conveniently  be  written  in  matrix  form  as  follows.  Define  r 
to  be  the  received-signal  vector  of  dimension  ,V, 


/  r(0)  \ 

/  r(U 

r  =  _  -  s  *  +  u  . 

\r(A~  1  )/ 

where  the  V-dimensional  vectors,  5r  and  w,  are  given  by 


/  MO)  \ 

/  «(0)  \ 

MO  1 

J 

a  ml 

u  = 

^  a  (  1  ) 

~  \)J  V"(  V-1  )/ 


-  8  - 


>  * 

,S,N 
s  v 

.■•/a 


7^v 

m 

.v>‘ 


*  vv 
.*v\ 

.  • 

A  •  « 

v.v. 


. 

•  •  • 

•  ./  »J 

v-.Vv 


;>r 


y.-$ 


Also,  define  S*  as  the  .V/Rx,V  rectangular  matrix  expressed  in  column-block  form  in  terms  of  /r 
separate  .Vx.V  matrices  according  to 


Vs-.  7 

where  S}  is  an  Vx.V  diagonal  matrix  containing  sample  values  of  the  complex  envelope  of  the 
transmitted  signal  sj(i). 


/ s T ( 0 , m *  j) 
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0  ■  ■ 
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'  0  s  T  ( 
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Further, 

define  the  reflectance 

vector  h 

of  dimension 

NI r  in  the  column- 

block 

form  of  /r 

vectors  according  to 
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where  each  Hi)  is  a  vector  of  dimension  V, 
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Using  (9)  and  these  definitions,  we  can  now  express  the  V-dimensional  signal  vector  vR  of  (14) 
and  (15)  as 


where  superscript  V  denotes  the  Hermitian-transpose  operation. 


In  terms  of  these  defined  matrices,  the  received  vector  has  zero  mean  and  covariance 


A  ,  -  E  r  r  '  »  E  s  Rs’R  *  E  uu. 


2E  T  S'  E  6b  *  5  +  \0l . 


Then,  since 


E>(t)6*(y):  = 

where  £(;)  is  the  Hermitian-symmetric  Toeplitz-matrix 


A(i)  = 


K  (0 ,  m  *  t)  K  *  (  1  ,  m  -*•  t )  •  •  A  *  ( .V  -  1  ,  m  *  i ) 

A  (  1  ,  rn  +  t )  K  (0  .  m  *  i) 


A  ( .V  -  I  ,  m  +  i) 


K  ( 0  .  m  +  t ) 


it  follows  from  (21)  that  the  covariance  A"r  of  r  is  given  by 

Kr“2FTS'KS*  ,V0/. 

where  K  is  the  block-diagonal  ,V/RX,V/R-dimensional  matrix  defined  by 


A  (  0 )  0  0 

0  A'  (  1  )  0 

0  0  0 


A  /,-  1 


The  \-ih  .diagonal  block  K(i)  of  K  is  the  covariance  matrix  of  the  reflectance  process  at  the  i -ih 
delav  bin. 


the  estimation  problem 


We  will  use  the  following  definition. 
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Definition:  Let  K  denote  the  set  of  all  ATrxAYr  block-diagonal  matrices  (25)  with  each  block 
K(i)  an  TVxTV  Hermitian-symmetric  Toeplitz  matrix  (23).  Let  fi  c  K  be  a  specified  convex  subset 
of  K.  Any  matrix  K  €  fi  is  termed  admissible.  A  variational  matrix  SK  e  K  is  called  an 
admissible  variation  of  K  for  a  fixed  K  €  f)  if  there  exists  an  a  >  0  such  that  K.  +  fiSK  e  fi  for 
all  fi  satisfying  0  <  fi  <  a. 

The  problem  is  to  form  an  admissible  estimate  of  the  covariance  matrix  K  of  (25)  given  the 
data  vector  r  of  (14).  The  radar  image  then  viewed  is  the  discrete  scattering  function,  an 
estimate  of  which  may  be  obtained  from  the  estimate  of  K  by  use  of  (13). 

In  the  definition  above,  the  constraint  that  K  be  in  fi  is  used  to  obtain  a  "reasonable"  setup 
of  the  problem.  Here,  we  assume  that  the  scattering  function  S(f,r)  in  (6)  is  only  nonzero  for 
frequencies  /  satisfying  |  /  |  <  /max  for  some  finite  upper  frequency  /max  and  for  all  delays 
This  is  equivalent  to  the  assumption  of  a  target  of  finite  cross  section  and  rotation  rate.  The 
discrete-time  scattering  function  S(v,i)  of  (13)  is  then  a  periodic  function  of  v  consisting  of  a 
sum  of  shifted  replicas  of  S(/,r)  scaled  in  amplitude  by  1  /At  and  in  frequency  by  A/,  where  A / 
is  the  time  between  samples  of  r(t).  The  replicas  of  5(/,r)  are  centered  at  every  integer  on  the 
v  scale.  In  order  to  guarantee  that  there  is  no  aliasing,  assume  that  the  sampling  rate,  //A/,  of 
r{t)  satisfies  the  Nyquist  condition,  //At  >  2/max.  Then,  S(vJ)  will  be  nonzero  between  -1/2 
and  +7/2  only  for  a  satisfying  1  v  |  <  i/max  =  The  output  of  our  algorithm  is  S(v,i) 

discretized  in  frequency  u.  For  a  resolution  having  at  least  1q r  (here,  the  subscript  CR  denotes 
cross  range)  samples  in  the  frequency  range  -^max  <  v  <  vmax,  a  total  of 


2  J7  /  max 

samples  of  u  between  -1/2  and  +7/2  are  required. 

The  model  for  the  incomplete  data  r  of  (14)  is  that  r  is  normally  distributed  with  zero 
mean  and  covariance  specified  in  (24).  Given  the  incomplete  data,  we  wish  to  estimate  the 
covariance  K  of  the  reflectance  process,  as  defined  in  (25).  To  do  this,  we  adopt  the 
maximum-likelihood  method  of  statistics,  which  selects  K  to  maximize  the  incomplete-data 
loglikelihood 

Lld{K)  =  -  -  In  (let  2EtS'  KS  *  A'/'  ’  -  -  r'f2FrS'  KS  *  V  ,/  ' '  r  . 

2  "2  (26) 
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where  the  maximization  is  subject  to  the  constraint  that  K  be  an  admissible  matrix. 

Lemma  1.  A  necessary  condition  for  an  admissible  K  in  the  interior  of  Cl  to  be  a  local 
maximum  of  Lt<i(K)  over  all  K  e  Cl  is 

Trfl  2£rS*  KS  +  ,V0/ rr'  -  2ErS'  KS-  N  0l){2E  TS'  K  S N0lV's'6KS}  =  0.  (27) 

for  all  admissible  variations  SK. 

The  proof  of  Lemma  1  in  the  Appendix  is  based  on  the  fact  that  the  necessary  condition 
for  an  admissible  K  to  maximize  L^(K)  is  that  for  all  admissible  variations  SK ,  there  holds 
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Lld{K  +  a6K)~  Lld(K)  ^ 
a 


<0. 


(28) 


We  call  (27)  the  trace  condition.  Burg,  Luenberger,  and  Wenger  [23]  have  studied  an 
equivalent  problem  of  Toeplitz-constrained  covariance-estimation  and  have  derived  the  trace 
condition  using  a  different  approach. 

If  ft  =  K,  there  are  A7R  unknowns  in  K.  Since  SK  e  K,  there  are  Wr  parameters  in  SK 
that  can  be  varied  for  this  case.  These  variations  in  the  trace  condition  generate  A7r  equations 
in  the  unknown  elements  of  K.  Thus,  in  principle,  the  trace  condition  produces  enough 
equations  to  determine  the  unconstrained  maximum-likelihood  estimate  K.  However,  the 
equations  are  complicated  due  to  the  inverse  matrices  appearing  in  (27),  so  it  does  not  appear  to 
be  feasible  to  determine  K  directly  from  the  trace  condition.  This  motivates  our  development 
of  the  iterative  approach  in  the  next  section.  A  sequence  of  estimates  that  increase  the 
likelihood  at  each  iteration  stage  is  identified,  and  we  demonstrate  that  limit  points  of  the 
iteration  satisfy  the  trace  condition. 

The  trace  condition  is  only  a  necessary  condition  which  the  estimate  K  must  satisfy.  For  it 
to  be  sufficient  as  well,  the  second  derivative  must  be  negative  along  all  admissible  variational 
directions  SK. 
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Lemma  2.  Sufficient  conditions  for  an  admissible  matrix  £  to  be  a  local  maximum  of  Lld(E) 
are  that,  first,  the  trace  condition  (27)  is  satisfied  for  all  admissible  variations  8K  and,  second, 
that  there  holds 

Tr[K’r'  S'  6KSK'r'{2ErS*  KS  +  N  al  -  2rr')K~l  S'  6KS}<0  (29) 

for  all  admissible  variations  5K 

The  proof  of  Lemma  2  is  in  the  Appendix.  Equation  (29)  is  just  the  second  derivative  of 


Lld(K )  in  the  direction  SK. 
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3.  Maximum- Likelihood  Imaging  for  the  Incomplete/Complete  Data  Model 

That  the  trace  condition  (27)  cannot  be  solved  directly  for  the  maximum-likelihood 
estimate  of  K  motivates  the  indirect  approach  we  now  take  of  embedding  the  imaging  problem 
in  a  larger,  seemingly  more  difficult  problem.  The  result  will  be  an  iterative  algorithm  which 
when  implemented  produces  a  sequence  of  admissible  matrices  K(°),  Kt1),  K(k),  ...  having  the 

property  that  the  corresponding  sequence  of  incomplete-data  loglikelihoods  Ljd[K(°)],  Z.id[K(i)], 
...,  Lid[Af(k)],  ...  is  nondecreasing  at  each  stage. 

Fuhrmann  and  Miller  [24]  have  recently  shown  that  maximum-likelihood  estimates  of 
Toeplitz-constrained  covariances  which  are  positive  definite  do  not  always  exist  when  given 
only  one  data  vector  r.  A  necessary  and  sufficient  condition  for  the  likelihood  function  to  be 
unbounded,  and  therefore  for  no  maximum-likelihood  estimate  to  exist,  is  that  there  be  a 
singular  Toeplitz  matrix  with  the  data  in  its  range  space.  For  our  imaging  problem,  this 
condition  is  that  there  exists  an  admissible  K  with 


singular  such  that 


2ETS' KS  +  N0I 


r  =  {2ETS’KS  +  N0l)a 


for  some  complex-valued  vector  a.  In  fact,  without  constraining  K  further  than  being  Toeplitz, 
a  sufficient  condition  that  a  singular  estimate  for  K  be  obtained  is  that  N0  =  0  and  there  exists  a 
singular  K  with  r  in  the  range  space  of  2EtS+KS.  The  argument  for  this  mirrors  that  of 
Fuhrmann  and  Miller  in  [24,  Theorem  1]  but  is  applied  to  the  complete  data  loglikelihood  given 
below  in  equation  (A7)  in  the  Appendix.  Furhmann  and  Miller  also  showed  that  even  if  the 
true  covariance  had  eigenvalues  bounded  from  above  and  below,  the  probability  that  there  exists 
a  singular  Toeplitz  matrix  with  the  data  in  its  range  can  be  very  close  to  one.  By  restricting  the 
search  to  Toeplitz  matrices  with  circulant  extensions,  they  were  able  to  show  that  the 
probability  a  singular  circulant  Toeplitz  matrix  has  the  data  in  its  range  space  is  zero.  Thus,  in 
order  for  maximum-likelihood  estimates  to  be  nonsingular  with  probability  one  for  all 
nonnegative  values  of  ,V0,  we  may  restrict  the  class  of  admissible  Toeplitz  matrices  to  be  those 
with  circulant  extensions  of  period  P,  where  P  is  equal  to  or  greater  than  the  number  ,V  of  data 
samples  available,  P  >  ,V.  What  we  envision  in  adopting  this  constraint  is  that  for  each  delay  /. 
the  .V  sample  values  of  the  reflectance  h(n.i),  n  =  0,  1 . V-/,  are  from  a  stationary  process 


.  -  V 


that  is  periodic  with  period  P,  where  P  could  be  some  large  but  finite  value;  a  lower  bound  on 
P  in  terms  of  a  desired  cross-range  resolution  is  discussed  above.  These  N  sample  values  of  the 
reflectance  enter  the  incomplete  data  r  according  to  (14)  and  (20).  By  using  the 
expectation-maximization  algorithm  of  Dempster,  Laird,  and  Rubin  [20],  we  shall  develop  a 
sequence  of  admissible  matrices  that  have  the  maximum-likelihood  estimate  of  K  subject  to  this 
circulant  extension  as  a  limit  point.  The  approach  parallels  that  of  Miller  and  Snyder  [21]  for 
estimating  the  power  spectrum  of  a  time-series  from  a  single  set  of  data.  An  important  benefit 
of  introducing  the  periodic  extension  and  using  the  expectation-maximization  algorithm  is  that 
estimates  of  both  the  scattering  function  and  the  reflectance  process  are  obtained  simultaneously 
and  can  be  readily  viewed  as  target  images  in  range  and  cross-range  coordinates:  thus,  the 
procedure  proposed  may  be  considered  to  be  natural  for  the  imaging  problem  because  both 
types  of  images  considered  separately  in  the  past  are  obtained  directly.  As  a  final  comment 
regarding  our  use  of  a  circulant  extension  for  K,  we  note  that  in  estimating  a  discretized  version 
of  the  target’s  scattering  function,  the  class  of  admissible  K  is  restricted  automatically  to  be 
those  with  circulant  extensions.  For  completeness,  we  also  include  in  the  Appendix  the 
equations  obtained  using  the  expectation-maximization  algorithm  for  estimating  general  Toeplitz 
matrices  when  the  assumption  of  a  circulant  extension  is  not  made. 

We  shall  introduce  a  modification  of  our  notation  to  indicate  that  the  ,V  samples  of  the 
reflectance  process  are  from  a  stationary  periodic-process  of  period  P.  Thus,  let  by(i)  denote 
the  Af-dimensional  vector  b(i)  of  (19).  We  now  think  of  by(i)  as  an  .V-dimensional  subvector  of 
the  /‘-dimensional  vector  bp(i)  formed  from  samples  of  the  reflectance  process  over  a  full 


period. 


MO  * 


6  ( 0  ,  m  *  i ) 
6(1  ,m  *  i) 


6(  V  -  1  .  m  -  () 


6  (  P  -  1  .  m  -  ; ) 


If  ly  is  the  Vx.Y  identity  matrix,  and  if  is  the  Px. V  matrix  defined  by 


-  is  - 


J  R 


then 


(32) 


b^O-J'/tbM. 

Also,  let  bfi  denote  the  A7-dimensional  vector  b  of  (18),  and  let  bp  be  the  ^/-dimensional  vector 
with  i-th  block  element  bp(i).  Then, 


b  v  \1  Rb P  , 

where  MR  is  the  PIx.NI  block  diagonal-matrix 


Mr- 


0  0 
Jr  0 


\0  0  0 


(33) 


Furthermore,  let  K^(i)  denote  the  NxN  Toeplitz  covariance-matrix  K(i)  of  AN(/)  defined  in 
(23),  and  let  Kp(i)  denote  the  PxP  circulant  covariance-matrix  of  bp(i).  Then,  the  Toeplitz 
matrix  Kn(/)  is  the  upper  left  block  of  the  circulant  matrix  Kp(i),  as  given  by 


Vv(0  J  R  ^  p(.  0  J  R  ' 

Lastly,  let  Kp  denote  the  PIRxPIR  block-diagonal  matrix  in  the  form  of  (25)  with  the  i-th 
diagonal  block  being  KP(i).  Then,  if  KN  denotes  the  NlRxNlR  matrix  K  of  (25),  there  holds 

A  y  =  M  p  A  p  \l  R  . 

Let  If'  denote  the  PxP  discrete  Fourier-transform  matrix  scaled  so  that  the  columns  are 
orthonormal. 


f 


.  o  o 

It  p  Up 


W  = 


\ 


0  * 
u  p  u  p 


o  p  - 1 
U  p  Up 
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it  , 


(P-  I  )* 


IV 


IP-  I  )(  P-  I  ) 


J 


(34) 


where  wP  =  exp(-j'2jr/P).  Also,  let  ft  p  be  the  PIRxPlR  block-diagonal  matrix 


16  - 


Then,  bp  can  be  represented  in  rotated  coordinates  according  to 


CL  p  Iv  p  t)  p 


f  a( 0)  \ 

/  a(D  \ 


'(36) 


\°  1  R~  0/ 

where  a(i)  *  Wbp(i).  The  assumption  that  bp(i)  originates  from  a  periodic  process  implies  that 
the  P/R-dimensional  vector  ap  is  normally  distributed  with  zero  mean  and  diagonalized 
covariance 


AP  =  E{aPa'P)  =  U'  Pk  rlv",,.  (37) 

We  will  denote  the  (p+ilp)-th  diagonal  element  of  Ap  by  a p2(i);  this  is  the  p-ih  diagonal  element 
of  the  Py-P  diagonal  matrix  E[a(/)a+(t)]. 

Let  S(v,i)  be  discretized  in  frequency  with  P  samples  taken  for  0  <  v  <1.  These  samples 
may  be  obtained  from  (13)  as 


5 


n  •  0 


K  F(n  ,i)e\ p 


2  nnp\ 

~J  P  )’ 


(38) 


for  p  =  0,  1 . P-1.  The  p-th  such  sample  is  just  the  p-th  entry  in  the  vector 

\  P  lv  K  P(i)e  , 

where  e  is  the  P-dimensional  unit  vector 


o  = 


(39) 


Substituting  (37)  into  (39),  we  get 


v  P  1/ K  P(i)e  =  v  P  A(i)  We  . 

But  P1/ We  is  a  P-dimensional  vector  of  ones  and,  therefore, 


(40) 


(41) 


which,  according  to  the  above  definition,  is  the  ( p+iln)-th  diagonal  element  of  the  diagonal 
matrix  Ap.  The  entries  of  the  diagonal  matrix  Ap  in  (37)  are  then  samples  of  the  scattering 
function. 

The  constraint  from  Section  2  that  the  scattering  function  S(v,i)  be  nonzero  only  for  M  < 
fmaxA*  =  ,  for  values  of  v  between  -1/2  and  + 1/2 ,  may  be  incorporated  at  this  point  in  the 

development.  Since  Ap  is  a  diagonal  matrix  of  samples  of  the  scattering  function,  we  restrict 
Ap  to  have  nonzero  values  only  in  its  top  left  and  bottom  right  corners.  More  precisely,  let  Iqr 
be  the  smallest  odd  integer  satisfying 

lc*  -  ■ 

Icr  is  the  number  of  cross-range  resolution  cells  implied  by  P  and  i/max.  Then,  let  Jcr  be  the 
IcrxP  matrix 


J 


CR 


/, 

0 


0 

0 


0 


where  Ii  is  an  [(/cR+7)/’M(/cR+f )/2\  identity  matrix,  and  /2  is  an  [(/cR-()/2]x[(fcR-^)/'?] 
identity  matrix.  Let  ,V/CR  be  the  IcrIrxPIr  matrix 


(Jcr  0  0  •  •  0  \ 

_  0  J  CR  0  0 

'•CR 

V  0  0  0  •  •  Jcj 

Define  Ep  to  be  the  diagonal  matrix 

P  =  M  CR  -4  P  M  CR  • 

The  diagonal  elements  of  EP  are  the  potentially  nonzero  diagonal  elements  of  .-Ip.  Recognizing 
that  some  elements  of  the  diagonal  matrix  .-tP  are  zero,  and  using  the  definition  of  AfCR,  we 
conclude  also  that 


-  18  - 


A  P  =  A/  CRZf  p  M  cr  . 

Then,  the  set  fl  referred  to  in  the  definition  in  Section  2  can  now  be  specified.  We  restrict 
consideration  to  those  covariance  matrices  generated  by  Ep  from  above,  so 

n  =  [K  =  M\W'pMcRL  pKi  CRW  PM  R). 

For  use  with  the  expectation-maximization  algorithm,  we  identify  the  complete  data  as 
(cp,w),  where  w  is  the  W-dimensional  noise  vector  defined  in  (15)  and  cp  is  defined  by 

C  P  ~  M  CRa  P  • 

Since  elements  of  np  corresponding  to  the  zero  diagonal  elements  of  A p  are  almost  surely  zero  , 
we  see  also  that 

a  p  =  M  crc  p  ■ 

Using  this  fact,  we  note  from  (14)  and  (20)  that  the  incomplete  data  r  can  be  obtained  from  the 
complete  data  according  to  the  mapping 

r-rcp  +  w,  (42) 

where  we  define  the  following  /r/crx^  matrix,  which  will  appear  throughout  the  development 
that  follows: 

r  =  v2£>  MCRWpMRS. 

The  loglikelihood  function  Lcd(Ep)  of  the  complete  data  as  a  function  of  Ep,  the  diagonal 
covariance-matrix  of  cP,  is  given  by 

,  !  '  l  ,  t  ■  ^  <  r-  '  \  1  -  v-  -  1 

L cd\£ p  j  p  j  j  ~  P  £ P  c p 

I*- 1  Ic«- 1  ,  !*-'  'cp-'  (43) 

=  £  in^co)-- £  I  :cp(t)  2o;2(i), 

l-O  p-0  z  c-0  p-o 

where  all  terms  that  are  not  a  function  of  Ep  have  been  suppressed,  and  where  cp(/')  is  the  p -th 
element  of  the  /CR-dimensional  vector  cp(/)  =  JcR*z(i)- 

The  expectation-maximization  algorithm  for  estimating  the  covariance  of  the  reflectance 
process  KP  from  the  incomplete  data  r  is  an  alternating  maximization  procedure  in  which  a 
sequence  of  estimates  of  Ep  having  increasing  likelihood  is  obtained  first.  If  Ep(*<)  denotes  the 
estimate  of  Ep  at  stage  k,  then  there  is  a  corresponding  element. 


Ky  =  w  F\i  CR£y'M  Cr\'p 

in  a  sequence  of  estimates  of  Kp  having  increasing  likelihood.  Likewise,  to  the  k-th  element 
Kp(k)  of  the  sequence  of  estimates  of  Kp,  there  is  a  corresponding  element, 

in  a  sequence  of  estimates  of  having  increasing  likelihood. 

Each  iteration  stage  of  the  expectation-maximization  algorithm  has  an  expectation  (E)  step 
and  a  maximization  (M)  step  that  must  be  performed  to  get  to  the  next  step.  The  E-step 
requires  evaluation  of  the  conditional  expectation  of  the  complete-data  loglikelihood  (43)  given 
the  incomplete  data  r  and  assuming  that  the  covariance  defining  the  complete  data  is  £p(k), 


Q[r,!r(/)]  =  E[Lcd(rf)|r,r(;)]. 


From  (43),  we  have  that 

/„-  I  /«"  1  I  >  CR~  1 

Q[r,!r">]-- I  I  m (o„(0)~I  I  E[|c,(i)|V.r<,*l]o;2<o-  ,45) 

i m  0  p  m  0  i m  0  p m  0 

The  M-step  yields  the  estimate  £P(k+1)  at  stage  k+1  as  the  choice  of  £P  that  maximizes  this 
conditional  expectation, 

I- p* * 1 5  =  a r g  max  [ Q ( ^ ^ 1 27 (/ } )] ,  (46) 

subject  to  the  constraint  that  the  maximizer  be  a  diagonal  covariance-matrix.  From  (45),  this 
maximization  yields  the  diagonal  matrix  Ep(k+i)  with  (p+ilp)-th  diagonal  element 


>p(0, 

Thus,  we  may  write  Ep(k+1)  as 


,(**!)_ 


E[  1  c p ( 0  | 2  r.r1;1]. 


r'p*'  n  =  E[  c  pc  r  .r'p0] . 

where  the  "d"  over  the  equal  sign  means  that  the  diagonal  terms  in  the  matrix  on  the  left  side 
equal  the  diagonal  terms  in  the  matrix  on  the  right  side  and  that  all  the  off  diagonal  elements 
on  the  left  side  are  zero. 


The  above  expression  (48)  appears  to  be  complicated  because  of  the  several  matrices  we 
have  defined,  but  it  produces  a  sequence  of  covariance  estimates  having  a  straightforward 
interpretation.  If  we  form  the  matrix  ATp(fc+i)  according  to 

a-  (P“ 1  >  =  w ;  \I  r * 1 5  .1/  CR  U ' ,, ,  (49) 

we  then  find  that 


A  p 


A*;*1'  (0) 


0  0 
A(t'n(l)  0 


V 


0 


0  0 


0  'N 

0 

a(;*,,:/r-i  J 


(50) 


a 


a 


■>' V 

#1 


where  Kp(k+1)(i)  is  a  PxP  circulant  matrix  interpreted  as  the  estimate  at  stage  k+1  of  the 
covariance  Kp(i)  of  the  P-periodic  reflectance  process  at  delay  m+i.  Miller  and  Snyder  [21] 
show  that  the  (n,m)-element  of  this  circulant  matrix  is  given  by 


P 


p- 1 

I 


p  -  0 


E[6(p.{)6*;<p 


+  m  -  n>  p 


a  ■(;>]. 


(51) 


where  <a> P  =  a  mod  P.  Equation  (51)  has  an  intuitively  appealing  form.  If  the  reflectivity 

process  b(n.i)  could  be  observed  for  all  instants  n  =  0.  1 .  P~1  in  a  period  and  for  each  i 

independently,  then  the  maximum-likelihood  estimate  of  the  covariance  Kp(i)  would  be  the 
arithmetic  average  of  the  lagged  products 

l'r'  ,  ,  <52) 

p  1  b(p,i)b*  <p  +  m-n>P, i  r 

'  p-0 


Equations  (50)  and  (51)  indicate  that  one  should  simply  substitute  the  conditional  mean  estimate 
of  an  unknown  lagged  product  into  this  expression  to  form  the  maximum-likelihood  estimate  of 
the  covariance  when  only  the  incomplete  data  are  known. 


estimating  EP  and  KP 

The  maximum-likelihood  estimate  of  Ep  is  a  limit  point  of  the  sequence  defined  in  (48). 


The  terms  on  the  right  side  of  this  equation  can  be  evaluated  as  follows.  Let  the 
conditional-mean  estimate  of  cp  in  terms  of  the  incomplete  data  r  be  defined  at  stage  k  by 


E[cP  r  ,  r  /  ] . 


Then,  (48)  can  be  rewritten  in  the  form 


rr  '  r  T  (*  ’  1  +  r (k  )  r 1  *  >* 

t-L  vc  />  c  f  c  p  j  '  L  p  \  C  p  C  p 


Now  examination  of  (42)  shows  that  forming  the  conditional-mean  estimate  (53)  of  c p  from  r  is 
a  standard  problem  in  linear  estimation-theory.  From  Tretter  [25,  Ch.  14],  for  example,  we 


find  that 


c(/>*,  =  r</)r[r*r(/)r  +  iv0/]*'r. 


Furthermore,  the  first  term  on  the  right  side  in  (54)  is  the  covariance  of  the  estimation  error 
when  cp  is  estimated  from  r.  Also  from  Tretter  [25,  Ch.  14],  we  have 

Fi  r  r  -  r  *  * '  '  *  r 

t[Cf  C  p  C  p  C  p  r  ,L  P  J 

(56) 

=  4‘)-r‘'lr[rr(;)rvv0/]'lrr';1. 

In  summary,  the  following  steps  are  performed  to  produce  a  sequence  Ep(°),  EpU), 

...,  Ep(*0,  ...  of  estimates  of  EP  for  which  the  corresponding  sequence  of  likelihoods  is 
nondecreasing: 

1.  set  A .  =  0,  select  a  starting  estimate  Ep(°); 

2.  calculate  the  estimate  of  cp  according  to  (55); 

3.  calculate  the  error  covariance  according  to  (56); 

4.  update  the  estimate  of  Ep  according  to  (54); 

5.  if  "last  iteration"  then  stop,  else  replace  k  by  k+1  and  go  to  2. 

The  starting  value  in  step  1  can  be  any  positive-definite,  diagonal  covariance-matrix  of 
dimension  IrIcr^r^cr- 

A  sequence  of  estimates  of  KP  having  increasing  likelihood  is  obtained  from  the  sequence 
of  estimates  of  Ep  according  to  |49) 


■/v'vfvf-.fvV.'-.'.v'v's 


V.V.V.V 


«*■  >  V ‘  > ">  v »  .*  . 


forming  the  scattering- function  image 

From  (4|),  the  diagonal  elements  of  the  /R/cRx/R/cR-dimensional  matrix  Ep  are  sample 
values  of  the  scattering  function,  with  the  scattering  function  samples  at  delay  m+i  given  by  the 
/cr  diagonal  elements  of  the  \-th  /cR*/cR-dimensional  diagonal  block  of  Ep.  We  may, 
therefore,  simply  regard  £p(k)  as  the  stage  k  estimate  of  the  scattering  function.  The  stage-k 
scattering-function  image  of  the  target  in  range  ( i  coordinate)  and  cross  range  (p  coordinate) 
can  be  displayed  as  follows.  Let  £p(k)(/)  denote  the  i-th  /CRx/CR-dimensional  diagonal  block  of 
Ep(k),  and  denote  the  p- th  diagonal  element  of  £p(k)(/)  by  s(p,i)  -  [crp2( / )](k)  for  p  =  0,  1,  ... 
1 cr- 1  •  Then,  s(0,i)  is  displayed  at  range  m+i  and  cross  range  corresponding  to  a  doppler  shift 
of  zero;  s(/,i)  and  s(Icr-IJ)  are  displayed  at  range  m+i  and  cross  range  corresponding  to  a 
doppler  shift  of  v  =  l/P  and  u  =  -1/P,  respectively;  s(2,i)  and  s(Icr-2,i)  are  displayed  at  range 
m+i  and  cross  range  corresponding  to  a  doppler  shift  of  v  =*  2/P  and  v  =  -2/P,  respectively;  and 
so  forth,  with  \ip,i)  and  sUCR-p,i)  displayed  at  range  m+t  and  cross  range  corresponding  to  a 
doppler  shift  of  ±p  P  for  p  =  1,2 .  ( Icr-D/2  when  /Cr  is  odd. 

forming  the  reflectance-process  image 

It  is  interesting  to  note  that  the  k -th  stage  conditional-mean  estimate  of  cp,  given  the 
measurements  r  and  assuming  that  the  second-order  statistics  of  reflectance  are  given  by  the 
k-;/i  stage  estimate  of  the  scattering  function,  is  used  to  form  the  estimate  of  EP  at  stage  A.W 
when  the  expectation-maximization  a'gorithm  is  used.  This  estimate  is  of  very  much  interest  in 
its  own  right  because,  from  (36)  and  its  definition,  the  / cr  elements  of  cpi  1 1  are  the  potentially 
nonzero  Fourier-transform  coefficients  of  the  reflectance  process  hpiu.  The  target’s  reflectance 
image  at  stage  k  is  formed  by  placing  these  elements  at  range  m+i  and  cross  range  in  the  same 
manner  as  described  above  for  the  scattering-function  image 

convergence  :  ■  -me \ 

There  are  some  important  properties  of  the  iteration  sequence  (48)  which  are  worth 
mentioning  f  :r  t.  each  step  is  in  an  improving  direction.  This  is  shown  by  writing  (54)  out  as 


w  here 
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and  where 

A-<*)  =  r-r<*>r+  yo/  (59) 

is  the  k -th  estimate  of  the  covariance  Kr  of  r.  Next,  the  trace  condition  (27)  which  the 
maximum-likelihood  estimate  must  satisfy  is  reexamined.  From  the  assumption  of  the 
P-periodicity  of  the  reflectance  process  and  the  matrix  definitions  given,  the  admissible 
variations  &K  must  be  of  the  form 

6  A'  =  \fRW'P\fCRbL  \i  CRWP\1  R.  (60) 

Here,  5E  is  a  diagonal  matrix  of  the  same  dimensions  as  E.  The  trace  condition  (27)  then 
becomes 

2  F  T  ' '  Tr  \K ;  V'  r '  L  P  r  *  ,V0  /  -  rr '  ,K ; 1  r '  bL  r]  =  0  .  (61 ) 

Using  the  fact  that  Tr(/lfl)  =  Tr(BA)  and  evaluating  this  trace  at  the  k-th  iterate,  we  see  that  the 
trace  on  the  left  side  of  (61)  is  equal  to 

2  ET  ■'  Tr  'elk)bZ'  .  (62) 

According  to  (57),  Ep(k)  is  changed  at  each  stage  by  adding  the  diagonal  elements  of 


to  Ep(k).  Define 


6r(*!  = 


as  these  diagonal  elements.  Then,  evaluating  the  trace  at  this  variation  gives 


This  shows  that  the  va.iation 


T  r  0 1  n  <5  JT  m  >  0 


6rm 


is  in  an  improving  direction.  Furthermore,  we  are  guaranteed  that  the  incomplete-data 
loglikelihood  is  nondecreasing  as  a  result  of  the  M-step  of  the  expectation-maximization 


algorithm  because  at  this  step,  the  conditional  expectation  of  the  complete-data  Ioglikelihood 
given  the  incomplete  data  and  the  last  iterate  for  £p  is  maximized  over  Up.  As  shown  in  [20] 
and  [21],  this  implies  that  the  incomplete-data  Ioglikelihood  is  nondecreasing. 

Lemma  3.  Assume  that  N0  >  0  and  EP(°)  is  positive  definite.  Then: 

(i),  £p(k)  is  positive  definite  for  all  k;  and  (ii),  all  stable  points  satisfy  the  trace  condition  (27) 
for  all  admissible  variations  (60). 

The  proof  of  the  first  part  of  Lemma  3  is  in  the  Appendix.  For  the  second  part,  since  the 
diagonal  elements  of  EP(k)  are  positive,  (65)  holds  with  equality  if  and  only  if  the  diagonal 
elements  of  0(k)  are  zero.  Notice  that  if  £p(k+1)  =  £p(k),  then  the  diagonal  elements  (64)  are 
zero.  This  implies  that  the  diagonal  elements  of  6(k)  are  zero  and,  hence,  that 

Tr(e(t)6r)  =  0  (66) 

for  all  diagonal  5E.  Thus,  all  stable  points  satisfy  the  trace  condition  (27)  for  admissible 
variations. 


computational  considerations 

The  computations  required  to  produce  radar  images  with  our  method  are  specified  by  (54), 
(55),  and  (56).  The  number  of  iterations  of  these  equations  that  are  required  to  produce  an 
image  near  the  convergence  point  is  presently  unknown.  Our  experience  in  using  an  iterative 
algorithm  to  produce  maximum-likelihood  images  for  emission  tomography  suggests  that  50  to 
100  iterations  may  be  necessary,  but  this  is  only  a  guess  that  will  not  be  verified  until  some 
experiments  are  completed.  Some  form  of  specialized  processor  to  accomplish  each  iteration 
stage  efficiently  will  probably  be  needed  to  produce  images  in  practically  useful  times.  One 
possible  approach  is  the  following.  The  matrix  product 

r  -  •i'2E~r  \t  CRW  F\l  KS 

is  required  at  each  iteration  stage  and  does  not  change.  This  /r/crx>'' -dimensional  matrix  can, 
therefore,  be  computed  once  off  line,  stored,  and  then  used  as  needed  during  on-line 
computations.  Then,  at  iteration  stage  k,  the  following  on-line  computations  can  be  performed: 


1.  compute  the  ,Vx.V-dimensional  matrix  A  defined  by  A  =  r^EpCOT  +  V0/; 

2.  compute  the  -dimensional  matrix  B  defined  by  B  =  £p(k)I\ 


3.  compute  BA~lr  and  the  diagonal  elements  of  £p(k)  -  BA  lB +. 

The  computations  in  3  can  be  accomplished  in  about  4A4/r/cr-2  time  steps  using  the  systolic 
array  described  by  Comon  and  Robert  [26]  augmented,  as  they  suggest,  by  one  row  to 
accomplish  the  postmultiplication  of  BA -1  by  r  and  by  /r/cr  rows  to  accomplish  the 
postmultiplication  by  B*.  The  matrix  multiplications  in  1  and  2  for  determining  A  and  B  can 
also  be  performed  rapidly  on  a  systolic  array.  More  study  of  implementation  approaches  is 
needed,  but  it  does  not  appear  that  the  computational  complexity  of  our  new  imaging  algorithm 
needs  to  be  a  limitation  to  its  practical  use. 

The  choice  of  N,  I r  and  Iq r  is  important  for  the  computations.  These  parameters  are 
selected  to  achieve  a  desired  range  and  cross-range  resolution  and  are,  therefore,  problem 
dependent,  but  the  same  considerations  used  with  other  approaches  to  radar  imaging  can  be  used 
in  selecting  them.  Choosing  the  product  /cZcr  to  be  on  the  order  of  N  will,  in  some  sense, 
make  the  imaging  problem  well  defined  because  the  number  of  unknown  parameters  ZcZcr  that 
need  to  be  estimated  to  form  the  image  is  then  comparable  to  the  number  of  measurements  N. 
On  the  other  hand,  the  choice  of  P  is  unique  to  our  approach.  As  stated,  we  need  P  >  A\  but 
no  upper  limit  is  given.  In  [24],  it  is  shown  that  as  P  increases  towards  infinity  so  does  the 
maximum  value  of  the  incomplete-data  loglikelihood  function,  with  probability  one.  Thus,  P 
cannot  be  made  arbitrarily  large  from  a  theoretical  standpoint.  Any  computation  involving  a 
matrix  with  one  dimension  equal  to  P  can  be  performed  off  line. 


4.  Conclusions 


The  expressions  we  have  obtained  for  forming  images  of  diffuse,  fluctuating  radar  targets 
are  based  on  the  model  stated  in  Section  2.  The  target  reflectance  is  assumed  to  introduce 
wide-sense-stationary  uncorrelated-scattering  (WSSUS)  of  the  transmitted  signal  with  no  glint  or 
specular  components  being  present.  The  reflectance  process  is  assumed  to  be  a  WSSUS  Gaussian 
process  with  unknown  second-order  statistics  given  by  a  delay-dependent  covariance  or 
scattering  function.  Echos  of  the  transmitted  signal  are  received  from  all  the  reflecting  patches 
that  make  up  the  target,  with  each  patch  introducing  some  propagation  delay,  doppler  shift,  and 
random  amplitude-scaling  into  the  signal  it  reflects.  The  superposition  of  the  echos  from  all  the 
patches  is  received  in  additive  noise.  Thus,  the  reflectance  process  is  only  observed  indirectly 
following  a  linear  superposition  and  in  additive  noise,  so  neither  the  reflectance  process  nor  its 
second-order  statistics  are  known.  Target  images  are  made  by  displaying  estimates  of  either  the 
reflectance  process  or  its  second-order  statistics  (scattering  function)  based  on  processing  the 
received  signal.  In  Section  2,  we  derived  the  trace  condition  which  the  maximum-likelihood 
estimate  of  the  covariance  of  the  reflectance  must  satisfy,  and  we  concluded  that  this  condition 
is  too  complicated  to  solve  explicitly  for  the  estimate.  This  motivated  the  introduction  in 
Section  3  of  the  incomplete-complete  data  model  and  the  use  of  the  expectation-maximization 
algorithm,  which  results  in  a  sequence  of  estimates  of  the  scattering  function  having  increasing 
likelihood.  A  corresponding  sequence  of  estimates  of  the  reflectance  process  is  also  obtained. 

There  are  a  number  of  issues  yet  to  be  resolved  for  the  approach  to  radar  imaging  which 
we  have  presented.  One  of  the  most  important  is  resolving  how  glint  and  specular  components 
in  the  return  echos  should  be  modeled  and  accommodated  in  the  formation  of  the  images.  The 
selection  of  transmitted  signals  to  produce  good  images  is  an  important  subject  about  which 
little  study  has  been  made.  The  quality  of  target  images  obtained  with  our  new  approach  is  not 
known  at  present;  to  study  this  issue,  we  are  presently  implementing  a  computer  simulation  so 
that  comparisons  to  alternative  processing  strategies  can  be  made.  The  equations  we  have 
developed  are  computationally  demanding,  so  special  processing  architectures  will  be  important 
to  make  their  use  practical. 


5.  Acknowledgment 

We  are  grateful  to  Dr.  Richard  E.  Blahut  of  Owego  NY  for  reading  a  draft  of  this 

paper  and  making  several  helpful  suggestions  for  improving  it  and  to  Dr.  Jeffrey  Shapiro  of 
Cambridge,  MA  for  bringing  Reference  [9]  to  our  attention. 


28 


6.  References 


1.  D.  R.  Wehner,  High  Resolution  Radar ,  Artech  House,  Dedham,  MA,  1987. 

2.  D.  L.  Mensa,  High-Resolution  Radar  Imaging,  Artech  House,  Dedham,  MA,  1984. 

3.  M.  Bernfeld,  "CHIRP  Doppler  Radar,"  Proceedings  of  the  IEEE,  Vol.  72,  pp.  540-541, 
April  1984. 

4.  D.  L.  Snyder,  H.  J.  Whitehouse,  J.  T.  Wohlschlaeger,  and  R.  Lewis,  "A  New  Approach  to 
Radar,  Sonar  Imaging,"  Proc.  SPIE  Conference  on  Advanced  Algorithms  and  Architectures, 
Vol.  696,  San  Diego,  CA,  August  1986. 

5.  J.  L.  Walker,  "Range-Doppler  Imaging  of  Rotating  Objects,"  IEEE  Transactions  on 

Aerospace  and  Electronic  Systems,  Vol.  AES- 16,  No.  1,  pp.  23-52,  January  1980. 

6.  J.  Shapiro,  B.  A.  Capron,  and  R.  C.  Harney,  "Imaging  and  Target  Detection  with  a 

Hetrodyne-Reception  Optical  Radar,"  Applied  Optics,  Voi.  20,  No.  19,  pp.  3292-3313, 
October  1981. 

7.  H.  L.  Van  Trees,  Estimation.  Detection,  and  Modulation  Theory.  Vol.  3,  John  Wiley  and 
Sons,  NY,  1971. 

8.  V.  S.  Frost,  J.  A.  Stiles,  K.  S.  Shanmugan,  and  J.  C.  Holtzman,  "A  Model  for  Radar  Images 

and  Its  Application  to  Adaptive  Digital  Filtering  of  Multiplicative  Noise,"  IEEE 

Transactions  on  Pattern  Analysis  and  Machine  Intelligence,  Vol.  PAMI-4,  No.  2,  pp. 
643-652,  March  1982. 

9.  N.  T.  Gaarder,  "Scattering  Function  Estimation,"  IEEE  Trans,  on  Information  Theory,  Vol. 
IT- 14,  No.  5,  pp.  684-693,  September  1968. 

10.  P.  E.  Green,  "Radar  Measurements  of  Target  Scattering  Properties,"  Radar  Astronomy,  V. 
V.  Evans  and  T.  Hagfors,  Eds.,  McGraw-Hill,  New  York,  Ch.  1,  pp.  1-78,  1968. 

11.  T.  Kailath,  "Measurements  in  Time-Variant  Communication  Channels,"  IRE  Trans,  on 
Information  Theory,  Vol.  IT-8,  pp.  829-8236,  September  1962. 

12.  R.  G.  Gallager,  "Characterization  and  Measurement  of  Time-and-Frequency  Spread 
Channels,"  M.I.T.  Lincoln  Laboratory  Tech.  Rpt.  352,  Lexington,  MA,  April  1964. 

13.  T.  Hagfors,  "Some  Properties  of  Radio  Waves  Reflected  From  the  Moon  and  Their  Relation 
to  the  Lunar  Surface,"  J.  Geophys.  Res.,  Vol.  66,  p.  777,  1961. 

14.  T.  Hagfors,  "Measurement  of  Properties  of  Spread  Channels  by  the  Two-Frequency  Method 
with  Application  to  Radar  Astronomy,"  M.I.T.  Lincoln  Laboratory  Tech.  Rpt.  372, 
Lexington,  MA,  January  1965. 

15.  R.  Price.  "Maximum-Likelihood  Estimation  of  the  Correlation  Function  of  a  Threshold 
Signal  and  Its  Application  to  the  Measurement  of  the  Target  Scattering  Function  in  Radar 
Astronomy,"  M.I.T.  Lincoln  Laboratory  Group  Rpt.  34-G-J,  Lexington,  MA,  May  1962. 

16.  M.  V.  Levin,  "Estimation  of  the  Second-Order  Statistics  of  Randomly  Time-Varying  Linear 
Systems,"  M.I.T.  Lincoln  Laboratory  Group  Rpt.  34-G-7,  Lexington,  MA,  November  1962. 

17.  L.  G.  Abraham,  et  al..  "Tropospheric-Scatter  Propagation  Tests  Using  a  RAKE  Receiver," 
1965  IEEE  Communications  Convention  Record,  pp.  685-690.  1965. 

18.  B.  M.  Sifford,  el  al.,  "HF  Time-  and  Frequency-Dispersion  Effects:  Experimental 
Validation  of  an  FSK  Error  Rate  Model,"  Stanford  Research  Institute  Tech.  Rpt.  4,  Menlo 
Park.  CA,  March  1965. 

19.  B.  Reiffen,  "On  the  Measurement  of  Atmospheric  Multipath  and  Doppler  Spread  by  Passive 
Means."  M.I.T.  Lincoln  Laboratory  Tech.  Note  1965-6  G-66,  March  1965. 

20.  A.  D.  Dempster,  N.  M.  Laird,  and  D.  B.  Rubin.  "Maximum  Likelihood  from  Incomplete 
Data  Via  the  EM  Algorithm,"  J.  of  the  Royal  Statistical  Society,  Vol.  B,  39,  pp.  1-37.  1977. 


s'  V  S' 


'Jf'  j-  ' 


21.  M.  I.  Miller  and  D.  L.  Snyder,  "The  Role  of  Likelihood  and  Entropy  in  Incomplete-Data 
Problems:  Applications  to  Estimating  Point-Process  Intensities  and  Toeplitz-Constrained 
Covariances,"  Proceedings  of  the  IEEE,  Vol.  75,  No.  7,  July  1987. 

22.  M.  Turmon  and  M.  I.  Miller,  "Performance  Evaluation  of  Maximum-Likelihood  Estimates 
of  Toeplitz  Convariances  Generated  with  the  Expectation-Maximization  Algorithm,"  1987 
Conference  on  Information  Sciences  and  Systems,  Johns  Hopkins  University,  Baltimore, 
MD,  pp.  25-27,  1987. 

23.  J.  P.  Burg,  D.  G.  Luenberger,  and  D.  L.  Wenger,  "Estimation  of  Structured  Covariance 
Matrices,"  Proc.  IEEE,  Vol.  70,  No.  9,  pp.  963-974,  September  1982. 

24.  D.  R.  Fuhrmann  and  M.  I.  Miller,  "On  the  Existence  of  Positive  Definite  Maximum-Likeli¬ 
hood  Estimates  of  Structured  Covariance  Matrices,"  IEEE  Transactions  on  Information 
Theory,  in  review. 

25.  S.  A.  Tretter,  Introduction  to  Discrete-Time  Signal  Processing ,  John  Wiley  and  Sons,  New 
York,  1976. 

26.  P.  Comon  and  Y.  Robert,  "A  Systolic  Array  for  Computing  5/1-1,"  IEEE  Trans,  on 
Acoustics,  Speech  and  Signal  Processing,  Vol.  ASSP-35,  No.  6,  pp.  717-723,  June  1987. 


7.  Appendix 

proof  of  Lemma  1 . 

From  the  definition  of  the  loglikelihood  function  in  (26),  we  have 


A  *  abK  - L 


i  A 


(A  1  ) 


- — r’f  K  *  a2E  TS'  6KS  '  -  A  1 1 V  -  \n  del  K  r  +  a2  E  7  S' b  K  S  del  A,1  . 

2a  ^  7  r  >  2a 

where  Kr  is  the  covariance  of  the  incomplete  data  r  as  given  in  (26).  Examining  the  first  term 
on  the  right,  we  have  that  it  equals 


r*A';lf:  l  +  a2ETS'6KSK'rlY'  -  l)r 

2 a  r  ^  T  J  1  (A2) 

=  —  r~K~rla2ETS'bKSK~r'l  +  a2ETS'bKSK~.[r 
2a  r  '  T 

=  l-r‘  k;'2EtS'  6KSK~r'  r  +  0(a). 


Examining  the  second  term  on  the  right  in  (Al),  we  have 


-  —  In  del  '  l +  a2ETS' 6KSK;'')  - - — ln(det(/  *  aB)) 

2a  '  r  2a 


=  -  — ln:l+aTr(B)  +  ...*a',det(5) 

2a  v  (A3) 

=  _~Tr(B)  +  0(a), 


where  B  is  defined  in  the  first  equality.  For  any  K  e  n  to  be  a  local  maximum,  a  small 
variation  in  K  in  an  admissible  direction  cannot  increase  Lld(K ),  or 


1 1 rn  -  Lul(K  *  a 6 K)~  L,d(K)  <  0 
a.o-a 


( A4) 


for  all  admissible  variations  5K.  If  K  is  an  interior  point,  then  -&K.  is  also  an  admissible 
variation  and  (A4)  becomes  an  equality.  Substituting  (A2)  and  (A3),  we  get 


r'k~'2ETS'6KSk~[r~Tr[2ETS'6KSK~])  =  0, 


which  is  the  trace  condition  (27). 


proof  of  Lemma  2 


Suppose  that  K  satisfies  (27)  and  (29)  for  all  admissible  variations  6K.  We  now  show  that 
(29)  is  simply  the  second  derivative  of  Lid(K)  in  the  direction  SK  by  taking  the  limit 

Tr( ,2EtS‘[K  +  <x6K)S  +  N 0l)'1  [rr*  -  N 0I -2EtS'{K  +  a6K)s) 

x(  2  E  T S '  K S  +  a.2  E  T S  ‘  6 K S  +  N 0 / )  '  S ' 6 K S 
-  K~l\rr'  -  i\  0{  -  2E  TS'KS)Kr'S6KS) 

=  Tr  K;'2ETS'6KSK;'i2ETS'  KS  +  N0I  -2rr’)K'r'2E  TS'6KS).  (A6) 

Thus,  the  conditions  of  Lemma  2  are  the  standard  sufficient  conditions  for  a  point  K  to  be  the 
local  maximum  of  Lid(K).  Equation  (27)  says  that  the  first  derivative  of  Lld(K)  is  zero  at  K. 
Equation  (29)  says  that  the  second  derivative  is  negative  definite  along  admissible  variations 
from  K.  A  necessary  condition  for  K  to  be  a  relative  maximum  is  that  this  last  expression 
evaluated  at  K  be  equal  to  or  less  than  zero  for  all  admissible  variations  6K.  Under  the 
assumptions  in  Section  4,  admissible  variations  are  given  by  (60).  Substituting  (60)  into  (A5) 
and  evaluating  for  all  diagonal  matrices  6£  gives  the  second-order  necessary  condition. 

estimating  a  general  Toe  pi  it:  matrix 

In  Section  4,  we  derived  a  sequence  of  estimates  for  a  covariance  matrix  subject  to  the 
constraint  that  the  estimates  must  be  circulant  Toeplitz  matrices.  For  completeness,  we  develop 
and  discuss  in  this  section  the  equations  for  estimating  a  covariance  matrix  subject  to  the 
weaker  constraint  that  the  estimates  be  general  Toeplitz  matrices.  Similar  equations  for  other 
constraints  on  the  Toeplitz  matrices  are  easily  obtained  by  mimicking  the  steps  in  the  main  body 
of  this  paper. 

Let  the  complete  data  be  (h.w),  and  let  b  be  normally  distributed  with  zero  mean  and 
covariance  K ,  as  given  in  (27).  The  complete-data  loglikelihood  is 


& 


* 


-  32  ■ 


Lcd  kP  =  --ln  det  kP  ~^bPKP  b  r . 


where  all  terms  that  are  not  a  function  of  K p  have  been  suppressed.  Maximizing  this  function 
gives  the  trace  condition 

Trf K'x{bb'  -  k)K’'6K)  =  0,  <a8) 

which  the  maximum-likelihood  estimate  K  must  satisfy.  Performing  the  E  and  M  steps  of  the 
expectation-maximization  algorithm  yields  the  following  iteration  sequence  for  the  elements 
K(n.i).  n  =  0,  1 .  N-l,  of  the  covariance  matrix  K(i)  defined  in  (23): 

J  N-n-  i  (A9) 

K{k*  [)(n  ,t)  =  -  -  E[  2_.  b(j  ,m  +  i)b*(j  +  n,m  +  i)\r  ,K<'k)]. 

A i  —  n  ;.0 


In  matrix  form, 

K{k’n  =  Kik)  +  2Er  Kik)S  K[k)~'(  rr'  -  2ETS'  A'(t)S  -  K  0 1)  K(rk)' '  S '  K(k) , 


where 


K[k)  =  2E  TS'  K(k)  S  +  N  0I , 


(A10) 


(All) 


If  this  iteration  converges  to  a  stable  point,  then  the  trace  condition  is  satisfied  at  that  point,  as 
may  be  shown  by  using  the  same  arguments  as  in  Section  4.  It  is  worth  restating  that  the  reason 
this  iteration  is  not  recommended  here  is  that  the  probability  that  the  iteration  sequence 
generates  a  singular  estimate  for  K  approaches  one  as  /V  gets  large.  By  restricting  consideration 
to  Toeplitz  matrices  with  circulant  extensions,  the  loglikelihood  function  is  bounded  with 
probability  one  for  finite  extensions  and  a  positive  definite  K  is  generated  with  probability  one, 
as  proven  by  Fuhrmann  and  Miller  [24], 


proof  of  Lemma  J.  part  (i) 

Assume  that  the  initial  guess  Ep(°)  for  Ep  is  positive  definite  and  that  .V0  >  0.  We  will  now 
show  that  if  Ep(k)  is  positive  definite,  then  so  is  Ep(k+i),  and  thus,  by  induction,  EpW  is 
positive  definite  for  all  k.  One  key  to  following  this  derivation  is  the  matrix  identity 


B(l  -  l/?)~'  =  (/  +  BA)'1  B. 


( A 1 2 ) 


This  identity  is  used  to  rewrite  (57)  as 


r(/*n  =  h  (/V  0r(Pk)  r  r'  r(pk)  +  ^lrik) 

+  L(Pk)  Trr'  r'  Z{k))H' 

=  n  0{h  L(Pk)  r){r'  r^k)  h') 
+(wr‘‘,rr)|r7'r‘t)H‘) 
+  NlHL{Pk)H' . 


where  we  have  defined  H  according  to 


h  =  (r(Pk)rr'  +  n0i)~[  . 

Clearly,  all  the  diagonal  elements  of  (A  13)  are  greater  than  or  equal  to  zero.  To  show 
are  strictly  positive,  we  look  at  the  last  term  and  get  that  the  i-th  diagonal  element  is 


1  A1  CX~  1 


NlHZ("H')u-Nl  Y.  W„(r(/))  (/T 


1-0 

'a'ca-I 


■  > 

'/« 


-nI  Y  !(^)j2(r(; 


1-0 


(*)) 

hr 


(A  1 3) 


(A  14) 
that  they 

(A  1 5) 


which  is  clearly  positive  when  N0  >  0  since  H  is  invertible  and  all  diagonal  elements  of  £p(>0 


4.2  Appendix  2.  Parameter  Values  Used  in  Confidence* Weighted  Simulation  Experiments 


This  appendix  contains  a  table  updating  parameter  values  given  in  the  first  six-month 
progress  report  of  this  contract.  These  are  the  values  being  used  in  computer-simulation  studies 
of  the  confidence-weighted  algorithm. 


We  assume  the  following  for  the  gaussian  envelope  waveforms: 


Quantity 


Symbol  Value 


Distance  to  target: 

Number  of  Bursts: 

Number  of  Pulses/Burst:  n 
Center  Frequency: 
Wavelength : 

Speed  of  Light:  C 


R 
N 
1 

fc  3  GHz 
lambda  0.1  m 
3xl0~8  m/sec 


12.5  nautical  miles  =  23150  m 
128 


Given 

Target  Angular  Rotation  Rate: 
Downrange  Resolution: 
Crossrange  Resolution: 

Ratio  of  major  to  minor 
axis  lengths  for  desired 
ambiguity  function: 


w 

rd 

rc 


The  following  are  uniquely  determined  for  the  burst  sequence: 


Target  Diameter: 

D 

64  *rd 

Round  Trip  Time: 

Tr 

2*R/C 

Delay  Resolution: 

rt 

2* rd/C 

Frequency  Resolution: 

rf 

2*w*rc/lambda 

Constant  1: 

al 

2*pi*e 

Constant  2: 

a2 

2*pi/e 

The  following  are  uniquely  determined 

for  a  given  burst,  burst  (i 

anale  major  axis  makes  with 

the  positive  delay  axis: 

angle 

( i~l )/n*pi 

Parameter  1 • 

Pi 

al*sin( angle)  2  ^  a2’cos(. 

Parameter  2: 

p2 

(a2  -  al) *sin( 2*angle) 

FWHM  Burst  Duration: 

T 

sqrt  ( pl*rt/r f )/( 2*pi ) 

FM  Sweep  Rate: 

b 

p2/ ( 8  *pi *T~2 ) 

Burst  Bandwidth: 

BW 

b*T 

Total  Illumination  Time: 

Tt 

sum  3  *T(  i  )  +■  (N-l)*Tr 

Total  Aspect  Change: 

W 

w*Tt 

If  we  take  rc  =  rd  =  0.6  m,  e 

=  1. 3e8 , 

w  =  8.88xl0~-3  radians/sec 

we  get  the  following  representative  values: 

Target  Diameter: 

D 

38. 4  m 

Round  Trip  Time: 

Tr 

lc4  usee 

Delay  Resolution: 

rt 

4  nsec 

Frequency  Resolution: 

rf 

0.1066  Hz 

Constant  1: 

al 

2 . 6e8  *pi 

Constant  2: 

a2 

1 . 5e-8  Tpi 

Minimum  FWHM  Burst  Duration: 

Tmin 

16.0  nsec 

Maximum  FWHM  Burst  Duration: 

Tmax 

2.07  sec 

Minimum  FM  Sweep  Rate: 

bin  in 

15.6  Hz, sec 

Maximum  FM  Sweep  Rate: 

bmax 

3.41  GHZ  sec 

Maximum  Burst  Bandwidth: 

BWmax 

346.  MHz 

Total  Illumination  Time: 

Tt 

506.  sec 

Total  Aspect  Change: 

W 

257.  degrees 

<!» 


nw.v 


